*****************
**** Table 4 ****
*****************

use "$data/paymentdata_ready" , clear

keep if t==1

log using "$output/table 4", replace

reg postratiow1 treatmentidentityd2 treatmentidentityd3 treatmentidentityd4 teamdummy* firststickerdummy* , cluster(dayteam)

test treatmentidentityd2 = 0
test treatmentidentityd3 = 0
test treatmentidentityd4 = 0

test treatmentidentityd2 = treatmentidentityd3
test treatmentidentityd2 = treatmentidentityd4
test treatmentidentityd3 = treatmentidentityd4


reg postratiow1 treatmentidentityd2 treatmentidentityd3 treatmentidentityd4 teamdummy* firststickerdummy* preratiow1 lastdebt_demeanedw1 prepaymentsmadedummy* , cluster(dayteam)

test treatmentidentityd2 = 0
test treatmentidentityd3 = 0
test treatmentidentityd4 = 0

test treatmentidentityd2 = treatmentidentityd3
test treatmentidentityd2 = treatmentidentityd4
test treatmentidentityd3 = treatmentidentityd4

capture program drop Ey_boot
program define Ey_boot, eclass
twopm postratiow1 treatmentidentityd2 treatmentidentityd3 treatmentidentityd4 ///
firststickerdummy* teamdummy* preratiow1 lastdebt_demeanedw1 prepaymentsmadedummy* /// ///
, firstpart(probit) secondpart(regress) vce(cluster dayteam)
margins, dydx(treatmentidentityd2 treatmentidentityd3 treatmentidentityd4) nose post
end
bootstrap _b, seed(340) reps(1000): Ey_boot

test treatmentidentityd2 = 0
test treatmentidentityd3 = 0
test treatmentidentityd4 = 0

test treatmentidentityd2 = treatmentidentityd3
test treatmentidentityd2 = treatmentidentityd4
test treatmentidentityd3 = treatmentidentityd4




**** **** This is for estimating the RI - P values in the bottom panel
quietly {
* Ensure complete cases
drop if treatmentall == .

* Generate treatment dummies
forval k = 0/3 {
    gen treat`k' = treatmentidentity == `k' if !missing(treatmentidentity)
}

* Count number of clusters per treatment group
forval k = 0/3 {
    preserve
    keep if treatmentidentity == `k'
    keep dayteam
    duplicates drop
    count
    scalar n_cluster`k' = r(N)
    restore
}

* Estimate original treatment effects 


* Reduced model
reg postratiow1 treat1 treat2 treat3 teamdummy* firststickerdummy*
margins, dydx(treat1 treat2 treat3) nose post
scalar b_reduced1 = _b[treat1]
scalar b_reduced2 = _b[treat2]
scalar b_reduced3 = _b[treat3]

* Full model
reg postratiow1 treat1 treat2 treat3 teamdummy* firststickerdummy* preratiow1 prepaymentsmadedummy* lastdebt_demeanedw1
margins, dydx(treat1 treat2 treat3) nose post
scalar b_treat1 = _b[treat1]
scalar b_treat2 = _b[treat2]
scalar b_treat3 = _b[treat3]

* Two-part model
twopm postratiow1 treat1 treat2 treat3 teamdummy* firststickerdummy* preratiow1 prepaymentsmadedummyB* lastdebt_demeanedw1, firstpart(probit) secondpart(regress) vce(cluster dayteam)
margins, dydx(treat1 treat2 treat3) nose post
scalar b_twopm1 = _b[treat1]
scalar b_twopm2 = _b[treat2]
scalar b_twopm3 = _b[treat3]

* Randomization inference 

set seed 18421
mat A = J(1000,3,.) // Reduced
mat B = J(1000,3,.) // Full
mat C = J(1000,3,.) // Two-part

quietly forvalues row = 1/1000 {
    preserve

    set seed 18421`row'

    tempfile clusterrand
    keep dayteam
    duplicates drop
    gen random_cluster = runiform()
    save `clusterrand'

    restore
    preserve
    merge m:1 dayteam using `clusterrand', nogen

    sort random_cluster
    egen n = group(random_cluster) if random_cluster != .

    replace treatmentidentity = 0 if n <= scalar(n_cluster0)
    replace treatmentidentity = 1 if n > scalar(n_cluster0) & n <= (scalar(n_cluster0) + scalar(n_cluster1))
    replace treatmentidentity = 2 if n > (scalar(n_cluster0) + scalar(n_cluster1)) & n <= (scalar(n_cluster0) + scalar(n_cluster1) + scalar(n_cluster2))
    replace treatmentidentity = 3 if n > (scalar(n_cluster0) + scalar(n_cluster1) + scalar(n_cluster2)) & n <= (scalar(n_cluster0) + scalar(n_cluster1) + scalar(n_cluster2) + scalar(n_cluster3))

    tab treatmentidentity, gen(treatmentidentitydummy)

    * Reduced model
    quietly reg postratiow1 treatmentidentitydummy2 treatmentidentitydummy3 treatmentidentitydummy4 teamdummy* firststickerdummy*
    margins, dydx(treatmentidentitydummy2 treatmentidentitydummy3 treatmentidentitydummy4) nose post
    mat A[`row',1] = _b[treatmentidentitydummy2]
    mat A[`row',2] = _b[treatmentidentitydummy3]
    mat A[`row',3] = _b[treatmentidentitydummy4]

    * Full model
    quietly reg postratiow1 treatmentidentitydummy2 treatmentidentitydummy3 treatmentidentitydummy4 teamdummy* firststickerdummy* preratiow1 prepaymentsmadedummy* lastdebt_demeanedw1
    margins, dydx(treatmentidentitydummy2 treatmentidentitydummy3 treatmentidentitydummy4) nose post
    mat B[`row',1] = _b[treatmentidentitydummy2]
    mat B[`row',2] = _b[treatmentidentitydummy3]
    mat B[`row',3] = _b[treatmentidentitydummy4]

    * Two-part model
    quietly twopm postratiow1 treatmentidentitydummy2 treatmentidentitydummy3 treatmentidentitydummy4 teamdummy* firststickerdummy* preratiow1 prepaymentsmadedummyB* lastdebt_demeanedw1, ///
        firstpart(probit) secondpart(regress) vce(cluster dayteam)
    margins, dydx(treatmentidentitydummy2 treatmentidentitydummy3 treatmentidentitydummy4) nose post
    mat C[`row',1] = _b[treatmentidentitydummy2]
    mat C[`row',2] = _b[treatmentidentitydummy3]
    mat C[`row',3] = _b[treatmentidentitydummy4]

    drop random_cluster n treatmentidentitydummy*
    restore
    di `row'
}

* Convert and combine matrices 

* Reduced model
clear
svmat A
rename A1 reduced_treat1
rename A2 reduced_treat2
rename A3 reduced_treat3
gen obs = _n
tempfile temp_reduced
save `temp_reduced'

* Full model
clear
svmat B
rename B1 full_treat1
rename B2 full_treat2
rename B3 full_treat3
gen obs = _n
merge 1:1 obs using `temp_reduced', nogen

* Two-part model
tempfile temp_merged
save `temp_merged'
clear
svmat C
rename C1 twopm_treat1
rename C2 twopm_treat2
rename C3 twopm_treat3
gen obs = _n
merge 1:1 obs using `temp_merged', nogen


* === Calculate p-values ===

* Reduced model
gen p_reduced1 = (abs(reduced_treat1) > abs(b_reduced1))
gen p_reduced2 = (abs(reduced_treat2) > abs(b_reduced2))
gen p_reduced3 = (abs(reduced_treat3) > abs(b_reduced3))
gen p_reduced4 = (abs(reduced_treat1 - reduced_treat2) > abs(b_reduced1 - b_reduced2))
gen p_reduced5 = (abs(reduced_treat1 - reduced_treat3) > abs(b_reduced1 - b_reduced3))
gen p_reduced6 = (abs(reduced_treat2 - reduced_treat3) > abs(b_reduced2 - b_reduced3))

* Full model
gen p1 = (abs(full_treat1) > abs(b_treat1))
gen p2 = (abs(full_treat2) > abs(b_treat2))
gen p3 = (abs(full_treat3) > abs(b_treat3))
gen p4 = (abs(full_treat1 - full_treat2) > abs(b_treat1 - b_treat2))
gen p5 = (abs(full_treat1 - full_treat3) > abs(b_treat1 - b_treat3))
gen p6 = (abs(full_treat2 - full_treat3) > abs(b_treat2 - b_treat3))


* Two-part model
gen p_twopm1 = (abs(twopm_treat1) > abs(b_twopm1))
gen p_twopm2 = (abs(twopm_treat2) > abs(b_twopm2))
gen p_twopm3 = (abs(twopm_treat3) > abs(b_twopm3))
gen p_twopm4 = (abs(twopm_treat1 - twopm_treat2) > abs(b_twopm1 - b_twopm2))
gen p_twopm5 = (abs(twopm_treat1 - twopm_treat3) > abs(b_twopm1 - b_twopm3))
gen p_twopm6 = (abs(twopm_treat2 - twopm_treat3) > abs(b_twopm2 - b_twopm3))
}
* === Summary of all p-values
sum p*


log close
